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Abstract 

The present work is concerned with modeling the non-linear response of fiber rein- 
forced polymer laminates. Recent experimental data suggests that the non-linearity 
is not only caused by matrix cracking but also by matrix plasticity due to shear 
stresses. To capture the effects of those two mechanisms, a model combining a 
plasticity formulation with continuum damage has been developed to simulate the 
non-linear response of laminates under plane stress states. The model is used to 
compare the predicted behavior of various laminate lay-ups to experimental data 
from the literature by looking at the degradation of axial modulus and Poisson’s 
ratio of the laminates. The influence of residual curing stresses and in-situ effect on 
the predicted response is also investigated. 

It is shown that predictions of the combined damage / plasticity model, in gen- 
eral, correlate well with the experimental data. The test data shows that there 
are two different mechanisms that can have opposite effects on the degradation of 
the laminate Poisson’s ratio which is captured correctly by the damage / plasticity 
model. Residual curing stresses are found to have a minor influence on the predicted 
response for the cases considered here. Some open questions remain regarding the 
prediction of damage onset. 

Key words: Fiber Reinforced Laminates, Polymer Matrix Composites, 
Computational Mechanics, Non-Linear Material Response, Continuum Damage, 
Plasticity, Puck Failure Criterion. 



Notation 


Indices: 

1, 2, 3 ... ply coordinates (1 - fiber, 2 - transverse in-plane, 3 - out-of-plane direction) 

I, n, t ... fracture plane coordinates (I - fiber, n - normal, t - transverse direction) 
x, y ... global coordinates (x - loading direction, y - transverse to load, in plane) 


Roman letters: 

C d 

£<init 

E d 

incl 

^init 

Ei 

/e 

f-Uc,ply 

^Th 

LT Ic,ply 


f-UIc,ply 

^Th 

' J IIc,ply 


Ga 

G incl /-find 

nl — ^nt 

I 

k 

n 

P\n Pi 2 

S 

s 

Sis 

t 

t piy 

Y 

Yl 


compliance tensor of damaged ply 
compliance tensor of initial (undamaged) ply 
elasticity tensor of damaged ply 

elasticity tensor of inclusions in Mori-Tanaka formulation 

elasticity tensor of initial (undamaged) ply 

Young’s modulus in z-direction 

inclusion aspect ratio normal to fracture plane 

factor of effort 

critical energy release rate for intra-ply cracking in mode I 
critical energy release rate for intra-ply cracking in mode I 
including residual stresses 

critical energy release rate for intra-ply cracking in mode II 

critical energy release rate for intra-ply cracking in mode II 

including residual stresses 

shear modulus for ij-shear deformation 

inclusion shear modulus in the fracture plane 

fourth order identity matrix 

parameter of shear plasticity law 

exponent of shear plasticity law 

slope parameters for Puck failure criterion 

Eshelby tensor 

nominal shear strength 

in-situ shear strength of a ply cluster in a laminate 
thickness of a cluster of equally oriented plies 
thickness of one single ply 
nominal transverse tensile strength 

in-situ transverse tensile strength of a ply cluster in a laminate 
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Greek letters: 


pi 


P 

lij 

uh 

I IS 

7, 

pi 

7n-0 

£ 

£ pl 

K 

Md 

pi 

^mp 

Pi 

M 12 

Pi 

^23 

2^21 

<?fp 

cr 
rr 

a ij 

Uni 


FPF 


@ nip 
e q 

$fp 

£ 


lay-up angle of off-axis plies 

engineering shear strain on plane i in direction j 

shear strain under in-plane simple shear at failure (<712 = S ls ) 

plastic shear strain component in the fracture plane in direction i 

magnitude of plastic shear strain in the fracture plane 

strain tensor 

plastic strain tensor 

normal strain component on plane i in direction i 
damage evolution parameter 

damage parameter for shear stiffness recovery under compression 
parameter for influence of normal stress on shear plasticity 
influence parameter under in-plane simple shear, a 12 
influence parameter under out-of-plane simple shear, 023 
minor in-plane Poisson’s ratio {v 2 \ = U12E2/E1) 
angle between t-coordinate and fracture plane shear stress vector 
traction vector of the fracture plane 
ply stress tensor 

ply stress tensor at failure (i.e. when ply failure criterion is fulfilled) 

stress component on plane i in direction j 

shear stress component on the fracture plane in direction i 

total shear stress on fracture plane (projection of traction vector) 

equivalent fracture plane stress 

fracture plane angle predicted by Puck failure criterion 

damage state variable 

damage state variable at saturation 
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1 INTRODUCTION 


The use of fiber reinforced polymer (FRP) composites is increasingly popular 
in industries where lightweight design is beneficial. In order to achieve further 
weight reductions without compromising the reliability of composite parts, ac- 
curate prediction of the material response is essential. Due to their complex, 
hierarchical micro-structure, the load response of laminated composites is in- 
fluenced by a number of physical mechanisms. Ultimate failure of laminates 
is typically caused by fiber failure or delamination. Prior to failure, plastic 
deformation and cracking of the weaker matrix constituent may lead to non- 
linearity and can influence other failure modes through load redistribution and 
by creating local stress concentrations. 

Although matrix failure normally does not lead to ultimate laminate failure 
directly, the modeling of non-linearities caused by the matrix is important 
for two reasons. On the one hand, accurate modeling of load redistribution is 
necessary with respect to the influence of the matrix response on other failure 
modes such as fiber failure or delamination. On the other hand, the service- 
ability of a structure may be determined by criteria other than strength that 
depend on the matrix response, for example, if a maximum allowable defor- 
mation requirement has to be met or if matrix cracks cannot be tolerated 
(e.g. in pressure vessels). Consequently, the mechanisms of matrix-dominated 
behavior of laminates need to be understood and their effects captured ap- 
propriately by laminate models in order to improve predictions of the load 
response and failure of laminates. 

The modeling of non-linearities caused by the matrix has been very much 
focused on continuum damage mechanics [12, 21] where it is assumed that 
an accumulation of brittle matrix cracks is responsible for the observed non- 
linearity. Rather than looking at discrete cracks, however, continuum damage 
models simulate the response of a cracked ply by modifying the elastic prop- 
erties of the homogenized ply depending on some internal state variables. A 
number of continuum damage models for stiffness degradation due to ma- 
trix cracking under plane stress states have been presented in the literature 
[1, 2, 4, 7, 11, 20, 24, 26, 28]. These models have proven to be successful in pre- 
dicting the non-linear response when the damaged plies experience primarily 
tensile stresses perpendicular to the fiber direction. Under shear-dominated 
loading, however, comparisons between model predictions and experimental 
data have been less satisfactory. Recent research suggests that the non-linear 
response under shear dominated ply loads cannot be attributed to brittle 
mechanisms alone [16, 30, 31]. 

To model the non-linear response of composite plies under shear-dominated 
loading, a plasticity model has recently been proposed [25] and combined with 
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a damage model developed previously [22, 24, 26]. It is implemented with an 
extended version of classical lamination theory (CLT) to provide for analy- 
sis of multi-axial laminates under plane stress states including thermal and 
moisture effects. The combined damage / plasticity model provides significant 
improvements over the original brittle damage model as shown by comparing 
predictions of the two models to experimental data [25]. In addition to provid- 
ing better correlation with experimental data, the combined model captures 
the non-linear response of uni-directional (UD) laminates as well as residual 
strains after unloading, and it is able to explain discrepancies between the 
shear response that can be observed when using different test methods. 

In the present paper, the combined damage / plasticity model is used to inves- 
tigate effects of in-situ strength and residual stresses on the non-linear load 
response of laminates. First, the formulation of the combined model is briefly 
reviewed. Next, the influence of non-linear shear behavior on the predicted 
in-situ strength following a method proposed by Camanho et al. [6] is dis- 
cussed. Predictions of the combined model are finally compared to two series 
of experimental tests by Varna et al. [31, 32]. 


2 COMBINED DAMAGE / PLASTICITY MODEL 


The combined damage / plasticity model assumes that damage occurs in the 
form of brittle matrix cracks that span the whole thickness of a ply and lead 
to a degradation of the homogenized ply stiffness. Damage can only develop 
in plies embedded in a multi-axial laminate because the first matrix crack in a 
UD laminate corresponds to ultimate failure. Consequently, any non-linearity 
prior to damage onset in embedded plies and all non-linearity in UD laminates 
is attributed to plasticity. The constitutive equation of the combined model 
that relates the ply stress tensor, cr, to the ply strain tensor, e, is given by 

a = E d (e - £ pl ) , (1) 

where e pl is the plastic strain tensor defined by the plasticity model and E d 
is the elasticity tensor of a damaged ply given by the damage model. Both of 
these tensors can contribute to the non-linear response. 

The damage and plasticity formulations used herein are based on the Puck 
failure hypothesis for matrix dominated failure in fiber reinforced composites 
[17, 20, 23]. According to Puck, fracture occurs in a plane that is parallel to 
the fiber orientation and defined by a fracture plane angle, 9f p , as depicted in 
Fig. 1, left. For plane stress states, the fracture plane is perpendicular to the 
laminate plane (0f p = 0) under combinations of transverse tensile stresses and 
in-plane shear or moderate transverse compression and in-plane shear. For 
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high transverse compression combined with shear, the fracture plane angle is 
non-zero and can be computed analytically. 

The Puck criterion for plane stress (Puck 2D) [17, 20, 23] is used in the damage 
model to predict the onset and evolution of damage under multi-axial stress 
states and to compute the fracture plane in which damage accumulates. In 
the plasticity model, the accumulation of plastic strain is also assumed to be 
associated with the fracture plane predicted by Puck 2D. This assumption is 
supported by recent experimental work showing shear bands in UD-laminates 
under uniaxial compression [3] . According to this study, the shear bands have 
the same orientation as ply cracks that develop when load is increased further, 
which suggests that the shear bands are precursors of ply cracks. 

The damage/plasticity model is implemented in a stand-alone code combined 
with CLT to provide a tool for the non-linear analysis of multi-axial lami- 
nates. An extended version of CLT is used to allow for the consideration of 
plastic strains as well as the strains caused by moisture and thermal loads 
(e.g. [5]). The formulation of the damage and plasticity models is explained 
in the following sections. 

2.1 Plasticity Formulation 

The plasticity law assumes that plastic strains are caused by shear bands 
with the same orientation as the fracture plane predicted by Puck 2D. The 
plastic shear strain in that plane, 7 ^, is related to the shear stress acting on 

the fracture plane, a n ^ = + a^ t , which is the projection of the traction 

vector, < 7 f p , onto the fracture plane (see Fig. 1, right). The plastic shear strain 
is assumed to have the form 



Fig. 1. Definition of fracture plane and corresponding coordinate system, l-n-t, with 
regard to the ply coordinate system, 1-2-3, by fracture plane angle, 8{ p (left); trac- 
tions on the fracture plane for 8f p ^ 0 (right). 
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with plasticity parameters, k and n, and an equivalent stress, defined as 

a nil> ~ \^nip\ T F m/> <J nn j (3) 

to account for the influence of normal stress on the non-linear shear behavior 
that is observed in experiments [17, 19]. The factor /jA is interpolated from 
the corresponding parameters for in-plane and out-of-plane shear, /4 and /a 23, 
respectively, as 

= ^i2sin 2 (V’) + ^4cos 2 (^) , (4) 

where if: is the angle between the directions of cr nt and a n ^, (cf. Fig. 1, right). 
The parameters /4 and /4 are considered to be material parameters that need 
to be derived from experimental data. In general, the two parameters are not 
the same due to the different effect of the ply micro-geometry in longitudinal 
and transverse direction resulting in a different influence of normal stresses 
on the longitudinal and transverse shear response. The parameter /i 2:i can be 
determined from stress-strain data of a uni-axial transverse compression test 
on a UD-laminate. The factor should be derived from experimental data 
of tests with varying stress ratio a 22/012- A detailed discussion of parameter 
identification for the plasticity model as well as a method for estimating 
and ^23 when the necessary experimental data is unavailable is given in [25]. 

Finally, splitting the plastic shear strain 7^ into its two components, 7^ and 
7nt , and transformation to ply coordinates results in a strain tensor given by 

pl = f(0, 0,0, 7?2>0 , 0) t for 

1(0,4,4,7i P 2,0,0) t for 


0fp = O 


Qf P ^ 0 


2.2 Damage Formulation 


The elasticity tensor of a damaged ply, E d , is predicted by a continuum dam- 
age model presented in [22, 24, 26]. In that model, a scalar damage state vari- 
able, £, is introduced as a measure for the amount of damage in a ply. This 
damage state variable can increase with load but can never decrease. The load 
acting on a ply is quantified by a factor of effort, / E , which is determined from 

o- FPF /e = <r , (6) 

where cr is a given ply stress state and cr tPF is the corresponding failure stress 
state determined from the Puck 2D failure criterion. The damage state variable 
is related to the factor of effort by a damage evolution law of the form 

£ j 0 for 

£ wat [ 1 - exp (~ (/e( 4 K 2 ) ~ 1)2 ) for 


/e < 

/e > TT7 


1 

1+K 


( 7 ) 
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with one damage evolution parameter k. The maximum amount of damage 
that can be reached in a ply is given by the damage state variable at satu- 
ration, £ sat . The general shape of the evolution law function is chosen based 
on experimental observations regarding the increasing crack density with load 
(e.g. [14, 20]). In terms of continuum damage modeling, the development of 
cracks is equivalent to an increase of the damage state variable, £. According 
to the evolution law in Eq. 7, damage starts to develop when /e = is 
fulfilled. The evolution parameter k, therefore, determines the damage onset 
load and controls how quickly damage progresses with an increase of load. 
For example, the evolution law converges to the step function for k — 0 such 
that damage onset occurs at /e = 1 and the final damage state £ sat is reached 
instantly. 

The effect of a given damage state on the elastic response of a ply is predicted 
by the Mori- Tanaka method [13]. By this approach, the elasticity tensor of 
a material containing aligned ellipsoidal inclusions is computed as a function 
of the inclusion aspect ratio, e n , and the elastic properties of the inclusions 
and the surrounding material. In the damage model, the Mori-Tanaka method 
is employed using penny-shaped inclusions that are aligned with the fracture 
plane predicted by Puck 2D. Note that these inclusions are not intended to 
represent actual cracks in the material, rather they are used to derive the 
anisotropic elasticity tensor of the damaged material in a thermodynamically 
consistent way. Based on the formulation of Tandon and Weng [29], the com- 
pliance tensor of the damaged material as a function of £ is given by 

c d = (V)” 1 = 

| J - £ [(E incl - E init ) : ( S - £(S - I)) + £ init ] _1 : [E incl - £ init ] J : C init , 

( 8 ) 


where E mc] denotes the elasticity tensor of the fictitious inclusions, E uut and 
C imt are the elasticity and compliance tensors of the initial (undamaged) ma- 
terial, respectively, S is the Eshelby tensor, and I is the 4 th order identity 
matrix. The elastic properties assigned to the inclusions depend on the stress 
state. If the normal stress on the corresponding fracture plane, a nn , is tensile 
such that cracks would be open, the inclusions become voids with zero stiffness 
(£ incl = 0). For compressive normal stress, the properties of the inclusions are 
defined to be the same as those of the initial (undamaged) ply material ex- 
cept for reduced shear moduli in the fracture plane, Gf and G'f , which are 
computed as 

Gif = Gii n t cl = /XDkn| , (9) 

where the factor /r D is a material parameter accounting for shear stiffness 
recovery due to friction at the crack faces. 



In order to fully define the damage model, the four parameters k, £ sat , e n , 
and (I i) need to be defined. Damage typically progresses very quickly with 
increasing load. Reasonable values for the evolution parameter are therefore 
in the range of K = 0.01 — 0.05, which means that damage onset in the model 
occurs at 99% — 95% of the nominal failure load. The saturation state is 
normally set to £ sat = 0.2 for lack of better information and because the Mori- 
Tanaka approach becomes increasingly inaccurate for £ > 0.2. It should be 
noted, however, that the choice of £ sat is not particularly relevant because final 
failure, e.g. due to fiber failure, is typically reached before crack saturation 
occurs. The inclusion aspect ratio normal to the fracture plane, e n , which 
enters into the computation of the Eshelby tensor, is chosen to be very small 
such as to resemble a crack- like geometry. As long as e n < 0.01, the exact 
choice of e n has little effect on the damage model predictions. Little is known 
about the correct choice of /q> A value in the range of = 10 — 15 has 
previously yielded good results, however, a conservative approach would be to 
assume no shear stiffness recovery between the crack faces and setting // D = 0. 


2.3 In- situ strength with non-linear shear behavior 


It has been found in experiments that the transverse tensile and shear strengths 
of a ply embedded in a multi-axial laminate are higher than those of a unidi- 
rectional (UD) material (e.g. [8, 15]). This effect is commonly termed ’in-situ’ 
effect. It increases with decreasing thickness of a ply (or the number of equally 
oriented plies clustered together) and also depends on the location of a ply in 
the laminate (inner or outer ply). In the damage model, the in-situ effect can 
be taken into account by employing in-situ strengths, Yf and S- ls , rather than 
UD laminate strengths in the Puck 2D criterion which is used in the model to 
determine damage onset under multi-axial stress states. 


For the current work, an analytical solution proposed by Camanho et al. [6] 
is used to compute in-situ strengths. The in-situ solution for thin embedded 
plies assumes an initial flaw in the form of a crack whose size is equal to the 
thickness of a ply or cluster of equally oriented plies, t. The in-situ strength 
is assumed to correspond to the uni-axial stress at which this initial crack 
would start to grow parallel to the fiber direction. The start of crack growth is 
determined based on fracture mechanics. For transverse tension, the approach 
gives the in-situ strength as 


with 



K 

V2l\ 


4 G 


Ic.ply 


7ra 0 A| 


22 


and cio = 



. . . inner ply 
. . . outer ply 


(10) 
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where G'i c ,pi y is the critical energy release rate for mode I intra-laminar crack 
propagation in the fiber direction, and E \ , E 2 , and u 2 i = V 12 E 2 /E 1 are the 
ply’s longitudinal and transverse Young’s moduli and minor Poisson’s ratio. 


The in-situ shear strength has to take the non-linear shear response into ac- 
count and is determined following Ref. [6] from 


2 



012(712) dyi2 


4GlIc ply , 

— with cio 

na 0 


tj 2 ... inner ply 

t ... outer ply 


( 11 ) 


where Gn CjP i y denotes the mode II critical energy release rate associated with 
intra-laminar crack propagation parallel to the fiber direction, and 7^ refers 
to the shear strain at the crack propagation load that corresponds to ai 2 — Si s . 
The non-linear relation for in-plane simple shear given by the plasticity model 
(prior to damage onset) is determined from Eq. 2 for = <712 as 


el , pi a i2 . 

712 — 7i2 + 712 — 7^ l~ 

Lti2 



(12) 


From Eq. 12, it follows that 


d7l2 = (c7 + ^(T)" ' < 13 > 

Using Eq. 13 in the integral on the left hand side of Eq. 11 leads to 

^ q2 , n ora+1 _ 2Gnc,pl y ,-..s 

rt/^i ^is ' ( I 1 \ 7 n ^is * 

The in-situ shear strength, S - ls , is given by the real positive root of Eq. 14. If the 
exponent n is an odd positive integer, a closed form solution can be obtained 
for Eq. 14 which, in that case, has exactly one real positive root. Hahn and Tsai 
[9] proposed to use a third order polynomial to approximate the non-linear 
shear response. It has been found, however, that a higher exponent typically 
yields a better description of the non-linear shear response [25] . It is therefore 
suggested here to choose an exponent of n = 5 or n = 7. 


3 COMPARISON TO EXPERIMENTAL DATA 


The combined damage / plasticity model is used to simulate the load response 
of glass fiber / epoxy laminates with varying lay-up under uniaxial tension. 
Predictions for the degradation of axial modulus and Poisson’s ratio as func- 
tions of axial strain are compared to experimental data by Varna et al. [10, 
31, 32], 
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Table 1 

Material data of the glass fiber / epoxy material [10, 31, 32] and parameters used in 
the damage / plasticity model. 


Elastic and thermal properties 


Ei 

e 2 

g I2 

^12 

AT 

Oi2 — 

[GPa] 

[GPa] 

[GPa] 


[K] 

[l/K] 

44.73 

12.8 

5.8 

0.3 

-120 

1 E-5 


Plasticity and damage parameters 
n k nil k e n £ sat 

7 147.1 MPa 0 0.05 0.001 0.2 


3. 1 Model Parameters for the Glass Fiber / Epoxy Material tested 


The material used in the experiments is a toughened glass fiber / epoxy system 
(material specifications are not given in the references). The determination of 
model parameters for the material system is discussed in this section. The 
ply properties and model parameters used in the analyses are summarized in 
Table 1. Elastic and thermal properties are taken from [31] with AT refer- 
ring to the assumed temperature change from a stress free state to ambient 
temperature and ct 2 — cty denoting the difference between the coefficients of 
thermal expansion in longitudinal and transverse directions. 

Parameters for the plasticity formulation (Sec. 2.1) are determined from the 
non-linear shear response derived from tensile tests on angle ply laminates, 
(±/? 4 ) s , with two different lay-up angles (3 = 27° and (3 — 40° [31]. The 
experimental data is shown in Fig. 2, left, including the analytical curve fit for 
n — 7 and k — 147.1 MPa. The two lay-up angles lead to different stress ratios 
of (T 22 /o'i 2 = —0.06 for (3 — 27° and 022/012 — 0.3 for (3 — 40°. Since there 
is little difference between the data from these two laminates, it is assumed 
that the shear response is independent of transverse normal stresses, and the 
influence parameter for in-plane shear, /i 22 ■ is set to zero. The parameter ^23 
is not relevant for the test cases shown here since the predicted fracture plane 
angle is always zero, resulting in -0 = 90° for plane stress states (see Eq. 4 and 
Fig. 1). 

For the damage model (Sec. 2.2), the damage evolution parameter is chosen as 
k, — 0.05, which leads to a quick progression of damage as is typically observed 
in experiments, and the inclusion aspect ratio is selected to be very small, 
e n = 0.001, to resemble crack like voids. The damage variable at saturation 
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is set to £ sat = 0.2, however, this parameter is found to have no influence on 
the results presented since saturation is not reached within the strain range 
considered. 

The degradation of elastic ply properties resulting from the chosen set of dam- 
age parameters is shown in Fig. 2, right, for the 90° ply of a (02/904) s laminate 
under uniaxial tension. Damage onset occurs at 95% of the predicted ply fail- 
ure load as a result of k, — 0.05. For a smaller value of k, damage onset in 
Fig. 2, right, would shift to slightly higher strains (at most, to e 22 = 0.6% 
for k, — 0). The different degradation of elastic properties in the model is 
controlled by the inclusion aspect ratio e n . For a spherical inclusion (e n = 1), 
all four degradation curves would coincide. For the thin, crack-like voids used 
in the damage model, there is almost no change of E\ and zq 2 , a pronounced 
reduction of E 2 , and less severe degradation of &'i 2 (Fig. 2, right). This char- 
acteristic is consistent with analytical solutions for the stiffness degradation 
of a cracked ply (e.g. [7]). 

Ply strength values and critical energy release rates for the glass fiber com- 
posite are not given in [10, 31, 32], The mode I critical energy release rate 
can be back-computed from data of the (0 2 /904) s laminate test. According to 
[10], cracking of the 90° layers for the (0 2 /904) s layup starts at approximately 
e xx = 0.6% (where x denotes the loading direction and y is the in-plane trans- 
verse direction). If curing stresses are disregarded, this strain state corresponds 
to a transverse stress in the 90° plies computed via CLT of a 22 = 76 MPa. As- 
suming that this value represents the in-situ transverse tensile strength, 
of an 8-ply cluster with f ply = 0.144 mm, the mode I critical energy release 
rate can be computed by inverting Eq. 10 with a 0 = 8t ply /2 = 0.576 mm as 
Gic iP i y = 0.4kJ/nr. This value lies on the upper end of the typical range of 




Yp Pl [%] Transverse strain, e 22 [%] 


Fig. 2. Non-linear shear response from experimental data of two different laminate 
tests by Varna et al. [31] and analytical curve fit for the plasticity formulation (left); 
degradation of elastic ply properties predicted by the damage model (right) . 
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Table 2 

Critical energy release rates and Puck 2D slope parameters used in the analyses. 

r< r< /^<Th s~iTh „,t 

LT Ic,ply ( - T IIc.,ply ^Ic^ply ^Ilc.ply P 12 P 12 

[kJ/m 2 ] [kJ/m 2 ] [kJ/m 2 ] [kJ/m 2 ] 

0.4 0.8 0.51 1.02 0.3 0.25 


energy release rates determined by standard tests for mode I fracture. If the 
same procedure is followed including curing stresses in the CLT equation by 
assuming a thermal load AT = — 120 K, a failure strain of e xx = 0.6% corre- 
sponds to transverse stresses 0-22 — 85 MPa which leads to G p . h p , y = 0.51 kJ/m 2 
following Ecp 10. 

Since there is no experimental data available regarding shear failure, the mode 
II critical energy release rate is estimated as Gn c , p i y = 2 Gi c , p i y which corre- 
sponds to the typical Gn C)P i y /Gi C)P i y ratio for ply fracture for several carbon 
fiber / epoxy composites given in [27]. The values of G'n C;P i y and Gn Pply lead 
to 8-ply in-situ shear strengths of S is — 72.2 MPa and S'?* 1 = 73.5 MPa, re- 
spectively, which are consistent with typical values for glass fiber / epoxy. It is 
clear that these values are only a rough estimate and a variation of Gn c , P iy will 
influence the predicted damage onset in plies that are loaded mainly by shear 
stresses. For the examples shown here, a higher value of Gn c , P iy would shift 
damage onset in off-axis plies to higher strains. However, the non-linearity 
under shear loading is dominated by plasticity rather than damage and the 
effect is therefore not very significant. 

Finally, the slope parameters for the Puck 2D criterion, p\ 2 and Pi 2 > are chosen 
as p \ 2 = 0.3 and p\ 2 — 0.25 which correspond to the values suggested by 
Puck for glass fiber / epoxy materials [18]. Strength data for fiber failure and 
compressive failure is not relevant for the test cases considered here. The 
parameters related to ply strength predictions are summarized in Table 2. 


3.2 Laminate tests (±/3/904) s 


The first series of tested laminates has a lay-up of (±/3/904) s with four dif- 
ferent angles f3 — 0°, 15°, 30°, 40°. Three analyses are performed for each of 
the laminates and compared to experimental data [32] as shown in Figs. 3- 
6. In a first analysis, the in-situ effect is not taken into account and 8-ply 
strengths without residual stresses are used for all layers (i.e. Y t = 76 MPa and 
S — 72.2 MPa). The corresponding curves in Figs. 3-6 are labeled as ’nom- 
inal’. The other two analyses use in-situ strengths computed from Eqs. 10 
and 11 based on the thickness of each ply cluster and its location (inner 
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or outer). The curves denoted by ‘in-situ’ are computed without residual 
stresses (Gi c ,pi y = 0.4kJ/m 2 , G'n C;P i y = 0.8kJ/m 2 ), while the analyses for 
curves ‘AT = — 120 K’ include the effect of residual stresses by a superim- 
posed thermal load and using G£ h ply = 0.51 kJ/m 2 and G^ ply = 1.02 kJ/m 2 . 
The strength values used in the three analyses are summarized in Table 3. 

The degradation in Figs. 3-6 is caused primarily by transverse cracking in the 
90° plies. In fact, Varna et al. [32] assumed that cracking occurs only in those 
layers. For the (02/904) s laminate (Fig. 3), all three analyses yield exactly 
the same result since below 2% axial strain the 0° plies do not develop any 
matrix cracking and the predicted response of 90° plies is the same for all 
three analyses. 


(0 2 /90 4 ) s , VamaOl 


(0 2 /90 4 ) s , VamaOl 




Fig. 3. Results for lay-up (02/904) s ; axial modulus (left) and laminate Poisson’s 
ratio (right) normalized by their initial value; experimental data from [32]. 


(pml5/90 4 ) s , VamaOl 


(pml5/90 4 ) s , VamaOl 




Fig. 4. Results for lay-up (±15/904) s ; axial modulus (left) and laminate Poisson’s 
ratio (right) normalized by their initial value; experimental data from [32]. 
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For the (±15/904) s laminate (Fig. 4), the ‘nominal’ analysis predicts damage 
onset in the angle plies at approximately e xx = 1.5% in addition to damage of 
the 90° plies. This leads to a kink in the degradation curve of Poisson’s ratio 
while the effect on Young’s modulus degradation is only minimal. Similar 
observations can be made for the ‘nominal’ predictions of the (±30/904) s and 
(±40/904 ) s laminates (Figs. 5 and 6, resp.). For these two laminates, however, 
the onset of damage in the angle plies already occurs at lower strains and with 
slightly more effect on the predicted modulus degradation. 


The analyses using in-situ strengths do not predict any cracking of the outer 
layers for the first two laminates within the strain interval considered. For 
laminates (±30/904) s and (±40/904) s , the Poisson’s ratio curves show two 
kinks which correspond to damage onset in the +/? and the —(3 plies, respec- 


(pm30/90 4 ) s , VamaOl 


(pm30/90 4 ) s , VarnaOl 




Fig. 5. Results for lay-up (±30/904) s ; axial modulus (left) and laminate Poisson’s 
ratio (right) normalized by their initial value; experimental data from [32]. 


(pm40/90 4 ) s , VamaOl 


(pm40/90 4 ) s , VamaOl 




Fig. 6. Results for lay-up (±40/904) s ; axial modulus (left) and laminate Poisson’s 
ratio (right) normalized by their initial value; experimental data from [32]. 
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Table 3 

Nominal and in-situ strength values used in the analyses. 


ply thickness 
ply location 

t = 
inner 

f ply 

outer 

in-situ 
t = 2 f ply t 

outer 

= 4f ply 
inner 

t = 8i ply 
inner 

nominal 

all 

Y t [MPa] 

213 

151 

107 

107 

76 

76 

S [MPa] 

101 

91 

82 

82 

73 

73 




in-situ, AT = — 

120K 



Y t [MPa] 

241 

170 

120 

120 

85 


S [MPa] 

104 

95 

85 

85 

76 



tively. The fact that +/3 and —j3 plies do not start to crack at the same time 
is a result of the different in-situ strengths of outer and inner plies according 
to Eqs. 10 and 11. For all four lay-ups, there is little difference between the 
predictions with and without residual stresses. 

It is interesting to note that the ‘nominal’ predictions, in general, show bet- 
ter correlation to experiments than the ones that include the in-situ effect. 
It should be kept in mind, however, that damage onset in the angle plies de- 
pends very much on shear strength, especially for high values of lay-up angle 
f3, and that the shear strength is an estimated value. Regarding the in-situ 
predictions, there are two additional aspects that lead to further uncertainty. 
First, there may be some interaction between damage evolution in neighboring 
layers due to local stress concentrations near cracks in adjacent plies that can- 
not be accounted for in a continuum damage approach. Second, the laminate 
becomes unsymmetrical when damage evolves differently in +/3 and —/3 plies. 
In that case, the loading of a test specimen is different from the assumed uni- 
axial loading. Taking all these aspects into consideration, it is difficult to draw 
any firm conclusions on the applicability of the in-situ strength predictions. 
However, judging from the model predictions, it is likely that cracking of the 
angle plies contributes to the measured degradation in Figs. 5 and 6. 


3.3 Laminate tests (0/ ± Pa/Oi / 2 ) s 


The second series of tests reported in [31] was performed on laminates with 
a stacking sequence (0/ ± /?4/0i/ 2 ) s and five angles f3 — 90°, 70°, 55°, 40°, 25°. 
Non-linearity in these laminates originates only from the /3-plies, which experi- 
ence varying stress ratios, 0 - 22 / 0 - 12 , depending on the angle (3. Since the experi- 


16 



mental data only provides laminate stresses and strains, ply stress states have 
to be computed via CLT and depend on the assumed constitutive response. 
The ply loading paths in au - 0-22 stress space computed from the laminate 
strains are shown in Fig. 7. Loading paths assuming linear elastic plies are 
depicted as thick solid arrows. For lay-up angles (3 = 90°, 70° and 55°, the 
arrowheads indicate the stress states corresponding to the laminate strain at 
which cracking initiated in the experiments [31]. For f3 = 40° and f3 — 25°, no 
cracking was observed during the tests, and the arrowheads indicate the start 
of non-linearity in the experiments. 

According to the combined model, stress states with a small ratio of 022 / 012 , 
i.e. here for laminates with angles f3 = 55°, 40° and 25°, lead to significant 
plastic strain prior to damage onset. The ply stress states during loading 
predicted by the plasticity part of the model (i.e. suppressing the onset of 
damage) are shown by dashed lines in Fig. 7 and deviate significantly from 
the stress ratios computed for linear elastic plies as the amount of plasticity 
increases. The curve for f3 — 55° is plotted up to e xx = 1.134% which represents 
the average damage onset strain in tests of that lay-up. For /? = 40° and 
/ 3 — 25°, the dashed curves terminate at = 2% which constitutes the 
strain range tested without onset of cracking. In other words, the ply stress 
states given by the three dashed lines represent strain states that did not 
cause any cracking in the three corresponding experiments. For f3 — 90° and 
f3 — 70°, the plasticity model does not predict any non-linearity and would 
therefore result in the same loading path as given by the black arrows. 



Fig. 7. Loading paths of /3-plies in (0/ ± /?4/0 1 / 2 )s laminates subjected to uniaxial 
tension and theoretical ply failure stresses using in-situ strengths and Puck 2D 
failure criterion. 
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Also shown in Fig. 7 is the Puck failure envelope (using 4-ply in-situ strengths) 
which determines damage onset in the damage / plasticity model. For the spe- 
cial case of (3 — 90, the 8-ply in-situ strength, = 76 MPa, has to be applied. 
Comparing the stress states computed from the laminate strains at the on- 
set of cracking to the Puck failure envelope reveals some discrepancies. In the 
(3 = 90° and /? = 70° laminates, damage seems to develop prematurely, i.e. the 
arrows in Fig. 7 do not reach the theoretical failure stress. On the other hand, 
the predictions of the plasticity model for (3 — 55° and f3 — 40° suggest that, 
in these tests, the stress states at failure (ends of dashed lines) exceed the fail- 
ure envelope. Since the value of Gn CiP i y (and therefore S\ s ) is only an estimate 
due to the lack of experimental data, the underprediction of damage onset 
for f3 = 40° merely indicates that the actual shear strength is higher than 
the estimated one. For (3 — 90°, the deviation may be acceptable considering 
the typical scatter of experimental data. The discrepancy for (3 — 70° and 
/ 3 — 55°, however, seems too high to be explained by scatter. Especially the 
fact that damage onset is overpredicted for (3 — 70° but severely underpre- 
dicted for (3 — 55° is surprising. Therefore, it has to be concluded that damage 
onset under multiaxial stress states is not yet completely understood and will 
require further investigation. 


In Figs. 8-12, the degradation of axial modulus and laminate Poisson’s ratio 
normalized by their respective initial values are shown for all five laminates. 
The experimental data points are taken from [31]. Note that the degradation 
curves for f3 — 90°, 70°, and 55° in [31] are given as a function of crack density, 
which can be converted to axial strain by the corresponding analytical ex- 
pressions also provided in [31]. Model predictions are performed using in-situ 
strengths with and without residual stresses. For each of the cases f3 — 70°, 55°, 
and 40°, a third analysis labeled ‘best fit’ is performed in which the strengths 
IT and S ls are adjusted such that damage onset matches the experimental 


(0/90 8 /0 1/2 )s» Varna99 


(0/90 8 /0 1/2 ) s , Varna99 



O 






Fig. 8. Results for lay-up (0/90g/0 1 / 2 ) s ; axial modulus (left) and laminate Poisson’s 
ratio (right) normalized by their initial value; experimental data from [31]. 
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data for each test case. This means that for those ’best fit’ analyses, the fail- 
ure envelope depicted in Fig. 7 is modified by changing ±7 and S is such that 
the envelope passes through the end point of the respective loading curve in 
Fig. 7 (i.e. the arrowhead for (3 — 70° and the end points of the dashed lines 
for (3 = 55°, and [3 = 40°). The purpose of these adjustments is to see whether 
the degradation behavior is captured correctly when the uncertainty of ply 
strengths is factored out. 

As can be seen in Figs. 8-12, the correlation between test results and model 
predictions is very good except for the discrepancy in damage onset discussed 
above. Similarly to the first series of tests in Section 3.2, there is only small 
difference between predictions with and without residual stresses. Also in anal- 
ogy to observations in the previous section, it is found that the different ef- 

(0/pm704/0 1/2 ) s , Vama99 

1.4 
1.2 
1 

0.8 
0.6 
0.4 
0.2 
0 

0 0.5 1 1.5 2 0 °- 5 1 1-5 2 

£ n [%] e xx f%] 

Fig. 9 . Results for lay-up ( 0 / ± 704/04/2)8; axial modulus (left) and laminate Pois- 
son’s ratio (right) normalized by their initial value; experimental data from [ 31 ]. 

(0/pm55 4 /0 1/2 ) s , Varna99 

1.4 
1.2 
1 

$ 0.8 
X 

> 0.6 
0.4 
0.2 
0 

0 0.5 1 1.5 2 0 0.5 1 1.5 2 

e xx 1%] ta 1%] 

Fig. 10 . Results for lay-up ( 0 / ±554/04/2)8; axial modulus (left) and laminate Pois- 
son’s ratio (right) normalized by their initial value; experimental data from [ 31 ]. 



(O/pmSSyOjQh, Varna99 




(0/pm70 4 /0]p) s , Varna99 
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fects contributing to non-linearity are much more apparent in the degradation 
curves of Poisson’s ratio than those of axial modulus. The non-linearity for 
/ 3 — 90° and /3 = 70° is caused by damage only and leads to a decrease in 
axial modulus as well as in Poisson’s ratio (Figs. 8 and 9). For /3 = 40° and 
f3 = 25°, on the other hand, non-linearity is primarily due to plasticity which 
reduces the axial modulus but increases the laminate Poisson’s ratio (Figs. 11 
and 12). In the (3 — 55° laminate, damage and plasticity both contribute to 
the non-linear response (Fig. 10). Consequently, the Poisson’s ratio increases 
at first as a result of plasticity, but when damage and plastic strains accumu- 
late simultaneously, the Poisson’s ratio stays approximately constant, i.e. the 
opposing effects of damage and plasticity on Poisson’s ratio cancel each other 
out. 


(0/pm40 4 /0 1/2 ) s , Vama99 (0/pm40 4 /0 1/2 ) s , Vama99 



Fig. 11. Results for lay-up (0/ ±404/0 1 / 2 ) s ; axial modulus (left) and laminate Pois- 
son’s ratio (right) normalized by their initial value; experimental data from [31]. 


(0/pm25 4 /0 1/2 ) s , Varna99 (0/pm25 4 /0 1/2 ) s , Vama99 



Fig. 12. Results for lay-up (0/ ±254/0 1 / 2 )s; axial modulus (left) and laminate Pois- 
son’s ratio (right) normalized by their initial value; experimental data from [31]. 
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In general, stiffness degradation related to damage seems to be slightly under- 
estimated by the model (see Figs. 8 and 9 as well as Section 3.2) while plastic 
strains tend to be slightly overpredicted. However, the overall correlation with 
experimental data regarding non-linearity is satisfactory and the different ef- 
fects of damage and plasticity on the non-linear response are captured very 
well by the combined damage / plasticity model. The main issue that requires 
further investigation is the discrepancy between model and experiments re- 
garding damage onset under multi-axial stress states. 


4 CONCLUSIONS 


A ply-level model for fiber reinforced composites is proposed that combines 
plasticity and continuum damage mechanics to predict the non-linear response 
of polymer composite laminates. The Puck criterion for plane stress states is 
used for predicting damage onset in the model. To account for the in-situ 
effect in thin embedded plies, an analytical fracture mechanics based solution 
is adopted to compute in-situ strengths as a function of ply thickness. 

The proposed model is used to predict the load response of various laminates 
under uniaxial tension, and results are compared to experimental data from 
two series of tests from the literature. Parameters for the plasticity formulation 
of the model are identified from test data independent from the test data used 
in the comparisons. The damage parameters that control stiffness degradation 
due to damage are chosen as typical values for glass fiber materials. 

The first series of tests consists of a thick 90° layer embedded in various angle- 
ply sublaminates. While most of the non-linearity in these laminates is due to 
damage in the 90° plies, it is demonstrated by the analyses that in some cases 
damage in the sublaminates is likely to contribute to the non-linearity. The 
second test series investigates non-linearity due to multi-axial ply stress states 
by using laminates consisting of angle-ply sublaminates embedded between 
0° layers. The main challenge in these tests is found to be the prediction 
of damage onset for various stress ratios a 22 / <7 12 ■ Apart from the unresolved 
problem of damage onset, the different effects of plasticity and damage on the 
non-linear response are captured very well by the proposed model. 

The influence of residual stresses on predictions for both test series is investi- 
gated. It is found that residual stresses have little effect on the results, which 
is partly due to the fact that ply strengths are determined from the onset 
of cracking in an embedded layer. Therefore, residual stresses are implicitly 
taken into account to a certain extent even when they are not modeled di- 
rectly. When ply strengths are determined from UD laminate tests, it is to be 
expected that residual stresses have a significant influence on damage onset. 
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An interesting observation made from the comparison to experiments is that 
the degradation of Poisson’s ratio can give valuable additional information. For 
example, the additional degradation due to successive damage onset in several 
layers predicted for test cases of the first test series, is much more apparent in 
the degradation of Poisson’s ratio than in that of axial modulus. Furthermore, 
in the second series of tests, the Poisson’s ratio degradation clearly shows that 
there are two different mechanisms responsible for the non-linear response. One 
mechanism leads to an increase of Poisson’s ratio while the other mechanism 
causes a decrease. In the present model the two mechanisms are interpreted 
as matrix plasticity and matrix cracking. Since both mechanisms result in a 
reduction of the axial modulus, the two mechanisms cannot be distinguished 
by looking at the degradation of axial modulus only. 
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